library(xlsx)
library(plyr)
library(dplyr)
library(tidyr)
library(ARTool)
library(ggplot2)
library(ggsignif)

setwd("H://OneDrive//研究//3-小实验TU//5-论文//TU对本外竞争的影响//数据与分析")
                                                                     #设置数据读取与输出文件夹

##################
###  数据输入  ###
##################
rm(list=ls())
grow <- read.xlsx("20220712.xlsx",sheetName = "grow", encoding = "UTF-8")  #生长数据
grow$W <- as.factor(grow$W)
grow$U <- as.factor(grow$U)
grow$ID <- as.factor(grow$ID)
grow$PT <- as.factor(grow$PT)

leaf <- read.xlsx("20220712.xlsx",sheetName = "leaf", encoding = "UTF-8")  #叶片数据
leaf$W <- as.factor(leaf$W)
leaf$U <- as.factor(leaf$U)
leaf$ID <- as.factor(leaf$ID)
leaf$PT <- as.factor(leaf$PT)

ps <- read.xlsx("20220712.xlsx",sheetName = "ps", encoding = "UTF-8")  #光合数据
ps$W <- as.factor(ps$W)
ps$U <- as.factor(ps$U)
ps$ID <- as.factor(ps$ID)
ps$PT <- as.factor(ps$PT)

cnp <- read.xlsx("20220712.xlsx",sheetName = "cnp", encoding = "UTF-8")  #元素数据
cnp$W <- as.factor(cnp$W)
cnp$U <- as.factor(cnp$U)
cnp$ID <- as.factor(cnp$ID)
cnp$PT <- as.factor(cnp$PT)


##################
###  定义函数  ###
##################
summary.yb <- function(data, group, value, na.rm=FALSE, conf.interval=.95, .drop=TRUE) {
  library(plyr)
  library(dplyr)
  library(ARTool)
  result <- list()
  for (i in c(1:length(value))) {
    result[[i]] <- list()
    length2 <- function (x, na.rm=T) {
      if (na.rm) sum(!is.na(x))
      else       length(x)
    }
    datac <- ddply(data, group, .drop=T,
                   .fun = function(xx, col) {
                     c(N    = length2 (xx[[col]], na.rm=T),
                       mean = mean    (xx[[col]], na.rm=T),
                       sd   = sd      (xx[[col]], na.rm=T),
                       min  = min     (xx[[col]], na.rm=T),
                       max  = max     (xx[[col]], na.rm=T)
                     )
                   },
                   value[i]
    )
    datac$se <- datac$sd / sqrt(datac$N) 
    datac$down <- datac$mean - datac$se
    datac$up <- datac$mean + datac$se
    result[[i]][[1]] <- datac
    
    m <- art(eval(parse(text = paste(paste(value[i],"~"), paste(group, collapse = "*")))), data=data)  
                                                                     #根据需要进行多因素分析
    result[[i]][[2]] <- anova(m)                                     #查看交互分析结果
    result[[i]][[3]] <- art.con(m, eval(parse(text = paste(paste("~ "), paste(group, collapse = "*")))))%>%
      summary() %>%
      mutate(sig = symnum(p.value, corr = FALSE, na = FALSE,
                          cutpoints = c(0, 0.001, 0.01, 0.05, 0.1, 1),
                          symbols = c("***", "**", "*", ".", " ")
                          )
             )                                                       #查看各单元之间交互作用的具体结果
  }
  return(result)
}

transpose.yb <- function(list, name) {
  New <- list()
  for (i in c(1:length(list))) {
    for (j in c(1:length(list[[i]]))) {
      dataf <- cbind(data.frame(Index = rep(name[i], length(list[[i]][[j]][,1]))), list[[i]][[j]])
      if (i == 1) {
        New[[j]] <- dataf
      } else {
        New[[j]] <- rbind(New[[j]], dataf)
      }
    }
  }
  return(New)
}

##################
###  数据分析  ###
##################
#二级数据计算
grow$TB <- grow$AB + grow$UB                                         #总生物量
grow$RS <- grow$UB / grow$AB                                         #根冠比
for (i in c(1:length(grow[,1]))) {
  grow$RY[i] <- grow$TB[i] / mean(grow$TB[grow$Treat == grow$Treat[i] & grow$PT == grow$PT[i] & grow$re == grow$re[i] & grow$ID == "Mono"])
}                                                                    #相对产量
cnp$CN <- cnp$C / cnp$N                                              #碳氮比
cnp$CP <- cnp$C / cnp$P                                              #碳磷比
cnp$NP <- cnp$N / cnp$P                                              #氮磷比
  
#统计描述与非参数方差分析
group <- c("W", "U", "ID", "PT")
result_grow <- summary.yb(grow, group, c("H", "RL", "BS", "TB", "RS", "RY"))
result_leaf <- summary.yb(leaf, group, c("leaf_a", "leaf_m"))
result_ps <- summary.yb(ps, group, c("A", "Ci", "gs", "E", "WUE", "chl", "Fv.Fm"))
result_cnp <- summary.yb(cnp, group, c("CN", "CP", "NP"))


##################
###  结果输出  ###
##################
#表格输出
Table <- c(transpose.yb(result_grow, c("H", "RL", "BS", "TB", "RS", "RY")),
           transpose.yb(result_leaf, c("leaf_a", "leaf_m")),
           transpose.yb(result_ps, c("A", "Ci", "gs", "E", "WUE", "chl", "Fv.Fm")),
           transpose.yb(result_cnp, c("CN", "CP", "NP")))
for (i in c(1,4,7,10)) {
  Table[[i]]$Treat <- paste(Table[[i]]$W, Table[[i]]$U, sep = "-")
  Table[[i]]$Treat[Table[[i]]$Treat == "0-0"] <- "CK"
  Table[[i]]$Treat[Table[[i]]$Treat == "1-0"] <- "T"
  Table[[i]]$Treat[Table[[i]]$Treat == "0-1"] <- "U"
  Table[[i]]$Treat[Table[[i]]$Treat == "1-1"] <- "T*U"
  Table[[i]]$Treat <- factor(Table[[i]]$Treat, levels = c("CK", "T", "U", "T*U"))
  Table[[i]]$ID <- factor(Table[[i]]$ID, levels = c("Mono", "Mix"))
}
sheetname <- c("grow_SD", "grow_ANOVA", "grow_MC", "leaf_SD", "leaf_ANOVA", "leaf_MC",
               "ps_SD", "ps_ANOVA", "ps_MC", "cnp_SD", "cnp_ANOVA", "cnp_MC")
for (i in 1:12) {
  if (i == 1) {write.xlsx(Table[[i]], "result.xlsx",  sheetName = sheetname[i])}
  else {write.xlsx(Table[[i]], "result.xlsx",  sheetName = sheetname[i], append = TRUE)}
} 

#绘图
Data_fig <- list()
Data_fig[[11]] <- Table[[1]]
Data_fig[[21]] <- Table[[4]]
Data_fig[[31]] <- Table[[7]]
Data_fig[[41]] <- Table[[10]]

Data_fig[[12]] <- separate(Table[[3]], contrast, 
                          into = c("T1","U1", "ID1", "PT1", "T2","U2", "ID2", "PT2"), remove = FALSE)
Data_fig[[22]] <- separate(Table[[6]], contrast, 
                          into = c("T1","U1", "ID1", "PT1", "T2","U2", "ID2", "PT2"), remove = FALSE)
Data_fig[[32]] <- separate(Table[[9]], contrast, 
                          into = c("T1","U1", "ID1", "PT1", "T2","U2", "ID2", "PT2"), remove = FALSE)
Data_fig[[42]] <- separate(Table[[12]], contrast, 
                          into = c("T1","U1", "ID1", "PT1", "T2","U2", "ID2", "PT2"), remove = FALSE)

for (i in c(12,22,32,42)) {
  Data_fig[[i]] <- Data_fig[[i]][(Data_fig[[i]]$T1 == Data_fig[[i]]$T2 & 
                                    Data_fig[[i]]$U1 == Data_fig[[i]]$U2 & 
                                    Data_fig[[i]]$p.value < .05) & 
                                   (Data_fig[[i]]$ID1 == Data_fig[[i]]$ID2 | 
                                      Data_fig[[i]]$PT1 == Data_fig[[i]]$PT2), ]
  Data_fig[[i]]$Treat <- paste(Data_fig[[i]]$T1, Data_fig[[i]]$U1, sep = "-")
  Data_fig[[i]]$Treat[Data_fig[[i]]$Treat == "0-0"] <- "CK"
  Data_fig[[i]]$Treat[Data_fig[[i]]$Treat == "1-0"] <- "T"
  Data_fig[[i]]$Treat[Data_fig[[i]]$Treat == "0-1"] <- "U"
  Data_fig[[i]]$Treat[Data_fig[[i]]$Treat == "1-1"] <- "T*U"
  Data_fig[[i]]$Treat <- factor(Data_fig[[i]]$Treat, levels = c("CK", "T", "U", "T*U"))
  for (j in c(1:length(Data_fig[[i]]$Index))) {
    if (Data_fig[[i]]$ID1[j] == Data_fig[[i]]$ID2[j]) {
          Data_fig[[i]]$y[j] <- max(Data_fig[[i-1]]$up[Data_fig[[i-1]]$Treat == Data_fig[[i]]$Treat[j] &
                                                   Data_fig[[i-1]]$ID == Data_fig[[i]]$ID1[j] &
                                                   Data_fig[[i-1]]$PT == Data_fig[[i]]$PT1[j] &
                                                   Data_fig[[i-1]]$Index == Data_fig[[i]]$Index[j]],
                                    Data_fig[[i-1]]$up[Data_fig[[i-1]]$Treat == Data_fig[[i]]$Treat[j] &
                                                   Data_fig[[i-1]]$ID == Data_fig[[i]]$ID2[j] &
                                                   Data_fig[[i-1]]$PT == Data_fig[[i]]$PT2[j] &
                                                   Data_fig[[i-1]]$Index == Data_fig[[i]]$Index[j]]) +
            .06 * (max(Data_fig[[i-1]]$up[Data_fig[[i-1]]$Index == Data_fig[[i]]$Index[j]]) -
                   min(Data_fig[[i-1]]$down[Data_fig[[i-1]]$Index == Data_fig[[i]]$Index[j]]))
      } else {
        Data_fig[[i]]$y[j] <- max(Data_fig[[i-1]]$up[Data_fig[[i-1]]$Treat == Data_fig[[i]]$Treat[j] &
                                                     Data_fig[[i-1]]$ID == Data_fig[[i]]$ID1[j] &
                                                     Data_fig[[i-1]]$PT == Data_fig[[i]]$PT1[j] &
                                                     Data_fig[[i-1]]$Index == Data_fig[[i]]$Index[j]],
                                  Data_fig[[i-1]]$up[Data_fig[[i-1]]$Treat == Data_fig[[i]]$Treat[j] &
                                                     Data_fig[[i-1]]$ID == Data_fig[[i]]$ID2[j] &
                                                     Data_fig[[i-1]]$PT == Data_fig[[i]]$PT2[j] &
                                                     Data_fig[[i-1]]$Index == Data_fig[[i]]$Index[j]]) +
          .12 * (max(Data_fig[[i-1]]$up[Data_fig[[i-1]]$Index == Data_fig[[i]]$Index[j]]) -
                 min(Data_fig[[i-1]]$down[Data_fig[[i-1]]$Index == Data_fig[[i]]$Index[j]]))
      }

  }
  Data_fig[[i]]$ID1[Data_fig[[i]]$ID1 == "Mono"] <- 1
  Data_fig[[i]]$ID2[Data_fig[[i]]$ID2 == "Mono"] <- 1
  Data_fig[[i]]$ID1[Data_fig[[i]]$ID1 == "Mix"] <- 2
  Data_fig[[i]]$ID2[Data_fig[[i]]$ID2 == "Mix"] <- 2
  Data_fig[[i]]$PT1[Data_fig[[i]]$PT1 == "N"] <- -.125
  Data_fig[[i]]$PT2[Data_fig[[i]]$PT2 == "N"] <- -.125
  Data_fig[[i]]$PT1[Data_fig[[i]]$PT1 == "I"] <- .125
  Data_fig[[i]]$PT2[Data_fig[[i]]$PT2 == "I"] <- .125
  Data_fig[[i]]$ID1 <- as.numeric(Data_fig[[i]]$ID1)
  Data_fig[[i]]$ID2 <- as.numeric(Data_fig[[i]]$ID2)
  Data_fig[[i]]$PT1 <- as.numeric(Data_fig[[i]]$PT1)
  Data_fig[[i]]$PT2 <- as.numeric(Data_fig[[i]]$PT2)
  Data_fig[[i]]$xmin <- pmin(Data_fig[[i]]$ID1 + Data_fig[[i]]$PT1, Data_fig[[i]]$ID2 + Data_fig[[i]]$PT2)
  Data_fig[[i]]$xmax <- pmax(Data_fig[[i]]$ID1 + Data_fig[[i]]$PT1, Data_fig[[i]]$ID2 + Data_fig[[i]]$PT2)
}
Fig <- list()
for (i in c(1:4)) {
  Fig[[i]] <- ggplot(Data_fig[[10*i+1]])+
    geom_line(aes(x = ID, y = mean, group = PT, color = PT), linetype =2, position = position_dodge(width = -0.5))+
    geom_pointrange(aes(x = ID, y = mean, ymin = down, ymax = up, color = PT), position = position_dodge(width = -0.5))+
    geom_signif(
      data = Data_fig[[10*i+2]],
      aes(xmin = xmin, xmax = xmax, annotations = "", y_position = y), 
      manual = TRUE, vjust = .5
    )+
    geom_text(data = Data_fig[[10*i+2]], aes(x = (xmin + xmax)/2, y = y, label = sig))+
    facet_grid(Index~Treat,scales= "free")+
    theme_bw()+
    scale_y_continuous(expand = expansion(mult = c(.05, .15)))
}


###存图
ggsave(file = "Fig1.svg",plot = Fig[[1]], width = 10, height = 16)
ggsave(file = "Fig2.svg",plot = Fig[[2]], width = 10, height = 4)
ggsave(file = "Fig3.svg",plot = Fig[[3]], width = 10, height = 14)
ggsave(file = "Fig4.svg",plot = Fig[[4]], width = 10, height = 6)

ggsave(file = "Fig1.png",plot = Fig[[1]], width = 10, height = 12)
ggsave(file = "Fig2.png",plot = Fig[[2]], width = 10, height = 4)
ggsave(file = "Fig3.png",plot = Fig[[3]], width = 10, height = 14)
ggsave(file = "Fig4.png",plot = Fig[[4]], width = 10, height = 6)
